Comprehensive analysis and accurate quantification of unintended large gene modifications induced by CRISPR-Cas9 gene editing

Most genome editing analyses to date are based on quantifying small insertions and deletions. Here, we show that CRISPR-Cas9 genome editing can induce large gene modifications, such as deletions, insertions, and complex local rearrangements in different primary cells and cell lines. We analyzed large deletion events in hematopoietic stem and progenitor cells (HSPCs) using different methods, including clonal genotyping, droplet digital polymerase chain reaction, single-molecule real-time sequencing with unique molecular identifier, and long-amplicon sequencing assay. Our results show that large deletions of up to several thousand bases occur with high frequencies at the Cas9 on-target cut sites on the HBB (11.7 to 35.4%), HBG (14.3%), and BCL11A (13.2%) genes in HSPCs and the PD-1 (15.2%) gene in T cells. Our findings have important implications to advancing genome editing technologies for treating human diseases, because unintended large gene modifications may persist, thus altering the biological functions and reducing the available therapeutic alleles.

native PAGE gel and HPLC showing HbA production in WT HUDEP2 and HbS production in S-HUDEP2 clones. S-HUDEP2 clones with exclusive HbS production. Hemoglobin AFSC control contains mix of HbA, HbF, HbS and HbC (bottom-up).

Figure S2
Figure S2. Design and validation of ddPCR assay for the quantification of the HBB copy number. The EvaGreen-based ddPCR HBB copy number assay consists of two separate PCR reactions: a primer pair targeting HBB on chr11 and a primer pair targeting reference gene on chr12. Alleles containing large deletions or chromosomal rearrangements that remove primer binding site/s cannot produce HBB signals. Therefore, we checked the copy number of HBB relative to a reference gene (CACNA1C). (A) An HBB forward primer binds 68bp upstream of the cut-site, and a reverse primer binds 100bp downstream of the cut-site. (B) ddPCR amplitude plot for 3 representative clones showing droplets containing signals from HBB amplification. (C) ddPCR amplitude plot for 3 representative clones showing droplets containing signals from CACNA1C amplification. (D) The copy number of HBB was normalized by the copy number of CACNA1C to figure out the number of alleles (0, 1 or 2) carrying large deletion (LD) in each clone. For clone 1, the absence of the HBB copy number indicates two LD alleles. For clone 5, 100% of HBB copy number indicates the absence of LD. For clone 8, a 50% HBB copy number indicates one LD and one small INDEL (SI).   The short-range PCR (S-R PCR) amplify the 300 bp region spanning the cutsite, and alleles containing LD would fail to generate PCR product. The L-R PCR amplify the 5.5 kb region spanning the cut-site, and alleles containing a LD within the long-range PCR primer binding sites generate smaller PCR products.    Fig. S4. For clone 1, the absence of S-R PCR band and two size-shifted L-R PCR bands suggest LD/LD genotype. For clones 2-4, one expected sized band and one downward shifted L-R band suggests SI/LD genotype. Although clone 8 does not have a size-shifted L-R PCR band, this clone has SI/LD genotype based on the ddPCR allelic drop-off assay.   Table S2). We previously showed that Cas9 cutting induced DNA double-strand breaks (DSB) in HBB could be repaired using the homologous sequences from the δ-globin gene (HBD) as an endogenous template, resulting in SCD mutation correction. We found that two clones had homozygous SCD mutation correction mediated by HBD gene conversion, and six clones failed to amplify S-R PCR products. As expected, no genotype with large deletion (LD) could be identified by S-R PCR. In contrast, the combination of three assays (S-R NGS, L-R PCR and ddPCR) gave different genotypes of the single-cell clones: out of the 46 clones identified as homozygous SI genotype by S-R NGS, only 4 clones were indeed homozygous SI while 42 clones carried LD. We found that the six clones that failed to amplify S-R PCR product all had LD/LD genotype, and the two clones with HBD conversion had LD. As shown in (A), 28% large deletion alleles (i.e., 56 alleles) occurred in 50% of clones (50 clones), which is consistent with the Hardy-Weinberg predictions (48%). Therefore, the use of S-R NGS significantly overestimated the percentage of SI alleles (97.8%) compared with that identified using the combination of three assays (71%). More significantly, our results showed that 50% of the single-cell clones had LD in at least one allele, which caused a significant reduction of HBB copy numbers in gene-edited S-HUDEP2 cells as quantified by ddPCR. (B) S-R NGS significantly overestimated the percentage of SI alleles; 97.8% SI and 2.1% HBD alleles were incorrectly identified by S-R NGS compared to 71% SI, 28% LD, and 1% HBD alleles identified by the combination of three essays. S-HUDEP2 treated with R-66S HiFi RNP and corrective ssODN was assayed as single cell clones using Illumina short-read sequencing (NGS), L-R PCR, and ddPCR copy number assays for clonal genotypes. (A) Representative agarose gel image showing the S-R PCR and L-R PCR products of 8 clones. Clone 1 contains a combination of an upward size-shifted and a downward size-shifted banding pattern, suggesting an approximately 100bp large insertion and a small INDEL (LI/SI genotype). Although the insertion allele was amplified by S-R PCR, the NGS analysis failed to detect this particular allele, providing a genotype of near 100% 26bp deletion. In order to capture this large insertion, the alignment threshold was reduced from >80% to >70%. NGS analysis under the adjusted alignment conditions indicated 35% of reads corresponded to a     The nine plasmid templates were linearized and pooled with specific molar ratios, with 80% Template 9 and 20% of Templates 1-8 combined. (B) The relative percentages of Templates 1-9 in the pooled plasmid standard were quantified by duplex probe-based ddPCR using barcodespecific primers and reference primers. The synthetic DNA library was then used as the standard for a 3-step L-R PCR to generate UMI-tagged and barcoded PCR3 products, sequenced using SMRT-seq to quantify the percentages of Templates 1-9. Based on the aligned CCS reads, Template 9 in PCR3 product was 54.38%, significantly decreased from 79.9% quantified by ddPCR in the original sample as standard, largely due to PCR errors. (C) Using the aligned CCS reads, LDs of the same start position and size were clustered together to identify unique LD patterns, and each large deletion pattern was mapped relative to the Cas9 cut-site. The CCS reads contain false-positive LDs different from that in Templates 1-8. (D) LDs identified using UMI consensus reads were mapped relative to the Cas9 cut-site. Only Templates 1-8 were identified, demonstrating SMRT-seq with UMI can accurately quantify LDs without false positives.   Note that 21 UMI consensus sequences have the LD of the same size (267bp) and start position, accounting for 0.6% of the total UMI consensus sequences.          To understand the types of large intergenic modification missed by HBG1-specific sequencing, we amplified and sequenced the 10 kbp region, including HBG1 and HBG2. We observed intergenic LD extending further upstream of the cut site on HBG2 and/or downstream of the cut site on HBG1, removing both HBG1 and HBG2.

Figure S24
A B

R-66S_2_T2
WT -1 -12 -8 -3   We delivered R-66S gRNA complexed with HiFi Cas9 and WT Cas9, respectively, into SCD HSPCs from Donor#1 and compared gene modification rates at the HBB on-target site as well as the known off-target site OT18. HiFi Cas9 and WT Cas9 treated samples showed similar LD rates (30.3% vs. 31.5%) and intermediate deletion rates (6.2 vs. 6.3%) quantified by SMRT-seq with UMI. WT Cas9 showed a higher rate of large insertion than HiFi Cas9 (2.2% vs. 1.6%). The LD rate at the OT18 in the WT Cas9 treated sample was 3.9%.   The raw sequencing data from Illumina MiSeq were demultiplexed by bcl2fastq from Illumina and merged using FLASH (52). Merged reads were aligned to reference genome hg19 using BWA-MEM (53), and the reads that were not spanning the cut site were filtered out with SAMtools (55). The split reads were identified using BEDtools (56) and further processed to break-point based variant calling, while the small INDEL patterns were generated by CRISPResso2 (30) using the unsplit reads.  Figure S11 were processed by LongAmp-seq. (A) LongAmp-seq measured 70% Template 9 (corresponding to unmodified HBB) while the frequency of Template 9 in the original DNA standard used for the library prep was 79.9%, showing underestimate of the unmodified allele and overestimate of the LD-containing alleles by LongAmp-seq, largely due to having more PCR duplicates of Templates 1-8 compared to that of Template 9. (B) LongAmp-seq identified LDs (≥200 bp) were clustered based on their sizes and locations to identify unique LDs. Instead of UMI cluster size-based filtering used in SMRT-seq, LD patterns with below 0.01% read number of total aligned reads were considered background noise and filtered out, which removed the false-positive LDs. Each unique LD was mapped relative to the Cas9 on-target cut-site. LongAmp-seq was able to identify the correct alleles (Templates 1-9) presented in the original DNA standard.  , SMRT-seq with UMI filtered and consolidated (SMRT UMI) and LongAmp-seq. LD rates were consistently highest by SMRT-seq without filtering for UMI, followed by LongAmp-seq and SMRT-seq with UMI. The percentage of LDs obtained using LongAmp-seq (quantified as the number of reads containing unique LDs divided by the total reads) was compared to the LD allele frequency quantified by SMRT-seq using UMI consensus reads and showed excellent correlation, although without UMI-based correction of PCR bias and error, LongAmp-seq gave slightly higher LD rates.      The R-66WT RNP with sickle ssODN were delivered into the K562 erythroleukemia cell line and harvested gDNAs at different time points over 3 days post-delivery, and analyzed by LongAmpseq. We compared the size distribution of LDs over time and found that the repair of longer LDs was slower than shorter LDs.

Figure S37
Figure S37. Nanopore long-read sequencing bioinformatics pipeline for detecting structural variants. FASTQ sequence files were processed into FASTQ using Guppy Basecaller. We first used NGMLR (57) to map all long reads reference human genome hg19., and the reads that mapped to the HBB region were analyzed for deletions and insertions calling. The reads that could not be mapped by NGMLR were further aligned by BWA-MEM (53) and filtered by SAMtools to include the chimeric reads carrying potential large deletions. The insertion profile was from NGMLR calling, and the large deletion profile included both NGMLR called reads, and BWA-MEM identified reads.